#!/usr/bin/perl -w
use strict;
use Getopt::Long;
use String::CRC;

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Usage

#88888888888888888888888888888888888888888888888888888888888888888888888888888888

my $usage = <<'USAGE';

USAGE:

    panda_sprot_to_fa
                         swissprot.dat
                         pubmed.list
                         [--every = 1000]

            swissprot.dat is the SWISSProt file you wish to process
            pubmed.list is a list of Accession codes linked to pubmed IDs
            The optional --every parameter sets the number of records between
                 each printed dot show the progress of the scan.

    This script converts SWISSProt data files into FASTA format, retaining
       the taxonomical NCBI ID.

USAGE


#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	Command line options

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
my $opt_every = 1000;
my $help;
GetOptions(	'every=i'	=> \$opt_every,
			'help'		=> \$help);

			
die $usage if ($help);
die "\nERROR:\n\tMissing Arguments:\n" .$usage
				unless (@ARGV >= 2);





my $progress_report = ($opt_every >= 0) ? $opt_every : 1000;


my $pubmedfile = pop @ARGV;
open (PUBMED, ">$pubmedfile");
binmode PUBMED;
binmode STDOUT;


	

	
#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	uniquify

#			removes duplicate entries in an array

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
sub uniquify(\@)
{
	my %hash;
	@hash{@{$_[0]}} = ();
	@{$_[0]} = keys %hash;	
	return @{$_[0]};
}	


#88888888888888888888888888888888888888888888888888888888888888888888888888888888

#	line_wrap

#			puts a newline every n letters
#			n.b. Do not use for constant parameters eg "$3"

#88888888888888888888888888888888888888888888888888888888888888888888888888888888
sub line_wrap($$)
{
	my @array;

	while (length $_[0])
	{
		push @array, substr $_[0], 0, $_[1], "";
	}
	$_[0]=join "\n", @array;
}
print STDERR "\t[" if ($progress_report);


my $reccnt = 0;
my $SQ = '';
my $DE = '';
my $OS;
my $AC = undef;
my @PubMed;
my $ID = undef;
my $TAX_LINE = '';


my $extra_count = 0;
my $line_cnt;

# read swissprot line by line
while (<>)
{
	$line_cnt++;
	chomp;

	# SP-TrEMBL ID
	if (/^ID/)
	{
		/^ID +([^ \t]+)/;
		$ID = $1;
	}

	# taxonomy
	elsif (/^OX/)
	{
		/^OX *(.*)/;
		$TAX_LINE .= $1;
	}

	# accession code (name)
	elsif (/^AC/gc)
	{
		/ *(\w+);/;
		$AC ||= $1;
	}

	# description
	elsif (/^DE/)
	{
		if (/^DE *(.*)/)
		{
			$DE .= ' ' if ($DE ne '');
			$DE .= $1;
		}	
	}

	# taxonomy names written out long (comma and "and" separated)
	elsif(/^OS/)
	{
		if (/^OS\s+(.+)/)
		{
			$OS .= $1;
		}	
		else
		{
			print STDERR "Non-conforming OS line # $line_cnt\n $_\n";
		}
		
	}

	# pubmed id
	elsif(/^RX/)
	{
		push(@PubMed, $1) if(/^RX.*PubMed=(\d+)/);
	}


	# sequence is always the last section of a swissprot record

	elsif (/^SQ/)
	{
		# check format validity
		die unless $SQ eq '';
		die "No accession code at line # $line_cnt" if !$AC;


		++$reccnt;
		print STDERR "." if ($progress_report && $reccnt % $progress_report == 0);

		# get sequence data
		while (<>)
		{
			last if (/^\/\//);
			chomp;
			tr/ //d;
			$SQ .= $_;
		}

			
		


		#get ncbi taxonomy IDs
			my @TAXID = ($TAX_LINE =~ /(\d+)/g);
			uniquify(@TAXID);
			die "ERROR!! no TAXID for $AC at line# $line_cnt" unless @TAXID;
			$extra_count += @TAXID - 1;


		#split os line into lots of square brackets
			#first remove all . & "and"s
			$OS =~ tr/.//d;
			$OS =~ s/and//g;
			#split on commas
			my @OS = split /,/, $OS;
			#add brackets
			$_ = "[".$_."]" foreach (@OS);
			#rejoin
			$OS = join " ", @OS;
		
					
		# create accession line consisting of permutations of accession codes and organisms
			my @acc_fragments;
			foreach my $taxa(@TAXID)
			{
				push(@acc_fragments, 'SP|'.$AC. '|NCBI_TAXID|'.$taxa.'|'.$ID.' '. $DE.$OS);
			}
			my $acc_line = join("\x01", @acc_fragments);

			print ">$acc_line\n", line_wrap($SQ, 60), "\n";
		
		# map accession codes to pubmed IDs.
			my ($crc_lo, $crc_hi) = crc ($SQ, 64);
			print PUBMED "$AC\t$crc_lo\t$crc_hi\t$_\n" foreach (@PubMed);

		# reset variables before next record
			@PubMed = ();
			$OS = '';
			$SQ = '';
			$DE = '';
			$AC = undef;
			$ID = undef;
			$TAX_LINE = '';
	}
	
}

print STDERR "]\n\t$reccnt sequences were extracted.\n";
print STDERR "\t", $extra_count, " extra fasta sequences generated because genes from\n".
		"\tmultiple organisms were mapped to the same swissProt record!\n"
				if ($extra_count);
